##################################################
## Author: Stephanie L. DeMora, Jennifer L. Merolla, Brian Newman, and Elizabeth J. Zechmeister
## Project: Republican Pushback on Patriotism-Linked COVID-19 Vaccine Messages: A Note on Moral Reframing
## Purpose: Replication File for Manuscript + SI Tables & Figures
## Date: January 7th 2025
##################################################


library(haven)
library(ggplot2)
library(broom)
library(knitr)
library(kableExtra)
library(dotwhisker)
library(sjPlot)
library(sjlabelled)
library(Rmisc)
library(stringr)
library(tidyverse)
library(questionr)
library(weights)
library(jtools)
library(margins)
library(interplot)
library(janitor)
library(nnet)
library(huxtable)

# Data:
df <- readRDS("Replication_Data.RDS")

# Set plot theme:
steph_theme <- list(
  hrbrthemes::theme_ipsum_rc(base_size = 12,grid= "",plot_title_size = 20,axis_text_size = 18),
  ggsci::scale_color_npg(),
  ggsci::scale_fill_npg(),
  theme(legend.position = "top", plot.caption = element_text(hjust = 0, face= "italic"),
  plot.title.position = "plot",
  plot.caption.position =  "plot"))


# ---------------------------------------------------------------------------------------------------------------
# ----------------------------   Manuscript Tables   ------------------------------------------------------------
# ---------------------------------------------------------------------------------------------------------------

# ---------------------------------------------------------------------------------------------------------------
# Table 1. Vaccine and Booster Intentions
# ---------------------------------------------------------------------------------------------------------------

mod1 <- lm(Q4A_Likely_to_Vaccinate100 ~ treatment, data = df, weights = WEIGHT)
mod2 <- lm(Q4B_Likely_to_SecondDose100 ~ treatment, data = df, weights = WEIGHT)
mod3 <- lm(Q4C_Likely_to_Boost100 ~ treatment,  data = df, weights = WEIGHT)
mod4 <- lm(Q4D_Likely_to_Boost_OfficialRec100 ~ treatment, data = df, weights = WEIGHT)

export_summs(mod1, mod2, mod3, mod4,
             model.names = c("Vaccinate","Second Dose","Boost", "Boost on Recommendation"),
             stars = c(`*` = 0.1, `**` = 0.05), error_format = "(SE = {std.error}, p = {p.value})",
             note = " * p < 0.1;  ** p < 0.05",
             coefs = c('(Intercept)' = '(Intercept)',
                       'Public Health Official' = 'treatmentPublic Health Official',
                       'Democratic Voter' = 'treatmentDemocratic Voter', 
                       'Republican Voter' = 'treatmentRepublican Voter'))

# ---------------------------------------------------------------------------------------------------------------
# Table 2. Attitudes Toward COVID-19 Vaccines
# ---------------------------------------------------------------------------------------------------------------

mod5 <- lm(Q5A_Benefits_GreaterThan_Risks ~ treatment, data = df, weights = WEIGHT)
mod6 <- lm(Q5B_Responsibility_to_Vaccinate ~ treatment, data = df, weights = WEIGHT)
mod10 <- lm(Q5C_Vaccines_Help_Economy ~ treatment, data = df, weights = WEIGHT)
mod7 <- lm(Q6_Family_Friend_Encourage ~ treatment, data = df, weights = WEIGHT)
mod8 <- lm(Q7_Support_Mandate_Caretakers ~ treatment, data = df, weights = WEIGHT)
mod9 <- lm(Q8_Support_Mandate_Universities ~ treatment, data = df, weights = WEIGHT)

export_summs(mod5, mod6, mod10, mod7, mod8, mod9,
             model.names = c("Benefits Outweigh","Responsibility","Help Economy", "Encourage", "Mandate Healthcare", "Mandate University"),
             stars = c(`*` = 0.1, `**` = 0.05), error_format = "(SE = {std.error}, p = {p.value})",
             note = " * p < 0.1;  ** p < 0.05",
             coefs = c('(Intercept)' = '(Intercept)',
                       'Public Health Official' = 'treatmentPublic Health Official',
                       'Democratic Voter' = 'treatmentDemocratic Voter', 
                       'Republican Voter' = 'treatmentRepublican Voter'))

# ---------------------------------------------------------------------------------------------------------------
# --------------------   SUPPLEMENTARY INFORMATION TABLES AND FIGURES   -----------------------------------------
# ---------------------------------------------------------------------------------------------------------------

# ---------------------------------------------------------------------------------------------------------------
# SI Table 1.  Balance Checks
# ---------------------------------------------------------------------------------------------------------------

regression <- multinom(treatment ~ GENDER + AGE4 + INCOME4 + EDUC5 + Ideology_num, data = df, na.rm = T, Hess = T, trace = F)
summary <- summary(regression)$coefficients/summary(regression)$standard.error
p <- (1 - pnorm(abs(summary), 0, 1)) * 2
p <- round(p, digits = 3)
p <- as.data.frame(p)
p$Comparison <- "Control"

#Compared to Republican Voter:
df %>%
  mutate(treatment = fct_relevel(df$treatment, "Republican Voter")) -> df_H2
regression <- multinom(treatment ~ GENDER + AGE4 + INCOME4 + EDUC5 + Ideology_num, data = df_H2, na.rm = T, Hess = T, trace = F)
summary <- summary(regression)$coefficients/summary(regression)$standard.errors
p2 <- (1 - pnorm(abs(summary), 0, 1)) * 2
p2 <- round(p2, digits = 3)
p2 <- as.data.frame(p2)
p2$Comparison <- "Republican Voter"

#Compared to Democratic Voter:
df %>%
  mutate(treatment = fct_relevel(df$treatment, "Democratic Voter")) -> df_H2
regression <- multinom(treatment ~ GENDER + AGE4 + INCOME4 + EDUC5 + Ideology_num, data = df_H2, na.rm = T, Hess = T, trace = F)
summary <- summary(regression)$coefficients/summary(regression)$standard.errors
p3 <- (1 - pnorm(abs(summary), 0, 1)) * 2
p3 <- round(p3, digits = 3)
p3 <- as.data.frame(p3)
p3$Comparison <- "Democratic Voter"

# Combine:
p_all <- bind_rows(p, p2, p3)
p_all %>%
  rownames_to_column("Comparison_2")

# ---------------------------------------------------------------------------------------------------------------
# SI Table 2. Descriptive Statistics by Experimental Condition
# ---------------------------------------------------------------------------------------------------------------

df %>%
  rename(Vaccinate = Q4A_Likely_to_Vaccinate100, 
         `2nd Dose` = Q4B_Likely_to_SecondDose100, 
         Boost = Q4C_Likely_to_Boost100, 
         `Boost on Recommendation` = Q4D_Likely_to_Boost_OfficialRec100,
         `Benefits Outweigh` = Q5A_Benefits_GreaterThan_Risks, 
         `Responsibility` = Q5B_Responsibility_to_Vaccinate,
         `Economy` = Q5C_Vaccines_Help_Economy,
         `Encourage` = Q6_Family_Friend_Encourage,
         `Mandate Healthcare` = Q7_Support_Mandate_Caretakers,
         `Mandate University` = Q8_Support_Mandate_Universities) %>%
  group_by(treatment) %>%
  rstatix::get_summary_stats(Vaccinate, 
                             `2nd Dose`, 
                             Boost, 
                             `Boost on Recommendation`,
                             `Benefits Outweigh`, 
                             Responsibility,
                             Economy,
                             Encourage,
                             `Mandate Healthcare`,
                             `Mandate University`,
                             type = "mean_sd") %>%
  arrange(variable)

df %>%
  tabyl(treatment)


# ---------------------------------------------------------------------------------------------------------------
# SI. Table 3. Vaccine Intentions controlling for intent to vaccinate
# ---------------------------------------------------------------------------------------------------------------

mod1 <- lm(Q4A_Likely_to_Vaccinate100 ~ treatment + Pre_Vax_Intentions, data = df, weights = WEIGHT)

huxreg('Likely Vaccinate' = mod1, 
       stars = c(`*` = 0.1, `**` = 0.05), error_format = "({std.error})", error_pos = "same",
       note = " * p < 0.1;  ** p < 0.05, standard errors in parentheses.",
       coefs = c('(Intercept)' = '(Intercept)',
                 'Democratic Voter' = 'treatmentDemocratic Voter', 
                 'Public Health Official' = 'treatmentPublic Health Official',
                 'Republican Voter' = 'treatmentRepublican Voter',
                 'Vaccination Intentions' = 'Pre_Vax_Intentions'),
       bold_signif = 0.1) 

# ---------------------------------------------------------------------------------------------------------------
# SI Table 4. Vaccine and Booster Intentions and Attitudes with Public Health Official set as baseline for comparison
# ---------------------------------------------------------------------------------------------------------------

df$treatment_2 <- relevel(df$treatment, ref = "Public Health Official")
mod1 <- lm(Q4A_Likely_to_Vaccinate100 ~ treatment_2, data = df, weights = WEIGHT)
mod2 <- lm(Q4B_Likely_to_SecondDose100 ~ treatment_2, data = df, weights = WEIGHT)
mod3 <- lm(Q4C_Likely_to_Boost100 ~ treatment_2,  data = df, weights = WEIGHT)
mod4 <- lm(Q4D_Likely_to_Boost_OfficialRec100 ~ treatment_2, data = df, weights = WEIGHT)
mod5 <- lm(Q5A_Benefits_GreaterThan_Risks ~ treatment_2, data = df, weights = WEIGHT)
mod6 <- lm(Q5B_Responsibility_to_Vaccinate ~ treatment_2, data = df, weights = WEIGHT)
mod7 <- lm(Q6_Family_Friend_Encourage ~ treatment_2, data = df, weights = WEIGHT)
mod8 <- lm(Q7_Support_Mandate_Caretakers ~ treatment_2, data = df, weights = WEIGHT)
mod9 <- lm(Q8_Support_Mandate_Universities ~ treatment_2, data = df, weights = WEIGHT)
mod10 <- lm(Q5C_Vaccines_Help_Economy ~ treatment_2, data = df, weights = WEIGHT)

huxreg('Likely Vaccinate' = mod1, 
       'Likely Second Dose' = mod2, 
       'Likely Boost' = mod3, 
       'Likely Boost on Rec' = mod4,
       'Benefits Outweigh' = mod5,
       'Responsibility' = mod6,
       'Encourage' = mod7,
       'Mandate Care' = mod8,
       'Mandate Uni' = mod9,
       'Economy' = mod10,
       stars = c(`*` = 0.1, `**` = 0.05), error_format = "({std.error})", error_pos = "same",
       note = " * p < 0.1;  ** p < 0.05, standard errors in parentheses.",
       coefs = c('(Intercept)' = '(Intercept)',
                 'Democratic Voter' = 'treatment_2Democratic Voter', 
                 'Republican Voter' = 'treatment_2Republican Voter',
                 'Control' = 'treatment_2Control'),
       bold_signif = 0.1) 

# ---------------------------------------------------------------------------------------------------------------
# SI Table 5. Attitudes Toward COVID-19 Vaccines Controlling for Vaccination Status
# ---------------------------------------------------------------------------------------------------------------

mod1 <- lm(Q4A_Likely_to_Vaccinate100 ~ treatment + DOV_VAXSTATUS_vaxdummy, data = df, weights = WEIGHT)
mod2 <- lm(Q4B_Likely_to_SecondDose100 ~ treatment + DOV_VAXSTATUS_vaxdummy, data = df, weights = WEIGHT)
mod3 <- lm(Q4C_Likely_to_Boost100 ~ treatment + DOV_VAXSTATUS_vaxdummy,  data = df, weights = WEIGHT)
mod4 <- lm(Q4D_Likely_to_Boost_OfficialRec100 ~ treatment + DOV_VAXSTATUS_vaxdummy, data = df, weights = WEIGHT)
mod5 <- lm(Q5A_Benefits_GreaterThan_Risks ~ treatment + DOV_VAXSTATUS_vaxdummy, data = df, weights = WEIGHT)
mod6 <- lm(Q5B_Responsibility_to_Vaccinate ~ treatment + DOV_VAXSTATUS_vaxdummy, data = df, weights = WEIGHT)
mod7 <- lm(Q6_Family_Friend_Encourage ~ treatment + DOV_VAXSTATUS_vaxdummy, data = df, weights = WEIGHT)
mod8 <- lm(Q7_Support_Mandate_Caretakers ~ treatment + DOV_VAXSTATUS_vaxdummy, data = df, weights = WEIGHT)
mod9 <- lm(Q8_Support_Mandate_Universities ~ treatment + DOV_VAXSTATUS_vaxdummy, data = df, weights = WEIGHT)
mod10 <- lm(Q5C_Vaccines_Help_Economy ~ treatment + DOV_VAXSTATUS_vaxdummy, data = df, weights = WEIGHT)


huxreg('Benefits Outweigh' = mod5,
       'Responsibility' = mod6,
       'Economy' = mod10,
       'Encourage' = mod7,
       'Mandate Care' = mod8,
       'Mandate Uni' = mod9,        
       stars = c(`+` = 0.2,`*` = 0.1, `**` = 0.05, `***` = 0.01),error_format = "({std.error})",
       coefs = c('(Intercept)' = '(Intercept)',
                 'Democratic Voter' = 'treatmentDemocratic Voter', 
                 'Public Health Official' = 'treatmentPublic Health Official',
                 'Republican Voter' = 'treatmentRepublican Voter',
                 'Vaccinated' = 'DOV_VAXSTATUS_vaxdummy'),
       bold_signif = 0.1) 


# ---------------------------------------------------------------------------------------------------------------
# SI Table 6. Vaccine and Booster Intentions and Attitudes with All Demographics as Controls
# ---------------------------------------------------------------------------------------------------------------

mod1 <- lm(Q4A_Likely_to_Vaccinate100 ~ treatment  + GENDER + AGE4 + INCOME4 + EDUC5 + Ideology_num, data = df, weights = WEIGHT)
mod2 <- lm(Q4B_Likely_to_SecondDose100 ~ treatment  + GENDER + AGE4 + INCOME4 + EDUC5 + Ideology_num, data = df, weights = WEIGHT)
mod3 <- lm(Q4C_Likely_to_Boost100 ~ treatment + GENDER + AGE4 + INCOME4 + EDUC5 + Ideology_num,  data = df, weights = WEIGHT)
mod4 <- lm(Q4D_Likely_to_Boost_OfficialRec100 ~ treatment  + GENDER + AGE4 + INCOME4 + EDUC5 + Ideology_num, data = df, weights = WEIGHT)
mod5 <- lm(Q5A_Benefits_GreaterThan_Risks ~ treatment  + GENDER + AGE4 + INCOME4 + EDUC5 + Ideology_num, data = df, weights = WEIGHT)
mod6 <- lm(Q5B_Responsibility_to_Vaccinate ~ treatment  + GENDER + AGE4 + INCOME4 + EDUC5 + Ideology_num, data = df, weights = WEIGHT)
mod7 <- lm(Q6_Family_Friend_Encourage ~ treatment  + GENDER + AGE4 + INCOME4 + EDUC5 + Ideology_num, data = df, weights = WEIGHT)
mod8 <- lm(Q7_Support_Mandate_Caretakers ~ treatment  + GENDER + AGE4 + INCOME4 + EDUC5 + Ideology_num, data = df, weights = WEIGHT)
mod9 <- lm(Q8_Support_Mandate_Universities ~ treatment  + GENDER + AGE4 + INCOME4 + EDUC5 + Ideology_num, data = df, weights = WEIGHT)
mod10 <- lm(Q5C_Vaccines_Help_Economy ~ treatment  + GENDER + AGE4 + INCOME4 + EDUC5 + Ideology_num, data = df, weights = WEIGHT)

huxreg('Likely Vaccinate' = mod1, 
       'Likely Second Dose' = mod2, 
       'Likely Boost' = mod3, 
       'Likely Boost on Rec' = mod4,
       'Benefits Outweigh' = mod5,
       'Responsibility' = mod6,
       'Encourage' = mod7,
       'Mandate Care' = mod8,
       'Mandate Uni' = mod9,
       'Economy' = mod10,
       #              stars = c(`*` = 0.1, `**` = 0.05), error_format = "({std.error})", error_pos = "same",
       #              note = " * p < 0.1;  ** p < 0.05, standard errors in parentheses.",
       coefs = c('(Intercept)' = '(Intercept)',
                 'Democratic Voter' = 'treatmentDemocratic Voter', 
                 'Public Health Official' = 'treatmentPublic Health Official',
                 'Republican Voter' = 'treatmentRepublican Voter',
                 'Gender' = 'GENDER',
                 'Age' = 'AGE4',
                 'Income' = 'INCOME4',
                 'Education' = 'EDUC5',
                 'Ideology' = 'Ideology_num'),
       bold_signif = 0.1) 


# ---------------------------------------------------------------------------------------------------------------
# SI Figure 1. Change in Vaccine Intentions and Attitudes Given Exposure to the Public Health Condition Relative to the Control Condition by Ideology (90% Confidence Intervals)
# ---------------------------------------------------------------------------------------------------------------

mod3 <- lm(Q4C_Likely_to_Boost100 ~ treatment * Ideology_num,  data = df, weights = WEIGHT)
mod7 <- lm(Q6_Family_Friend_Encourage ~ treatment * Ideology_num, data = df, weights = WEIGHT)
mod8 <- lm(Q7_Support_Mandate_Caretakers ~ treatment * Ideology_num, data = df, weights = WEIGHT)
mod9 <- lm(Q8_Support_Mandate_Universities ~ treatment * Ideology_num, data = df, weights = WEIGHT)

plot1 <- interplot(m = mod3, var1 = "treatment", var2 = "Ideology_num",facet_labs = c("Boost"),
                   point = T,ci = .90) +
  steph_theme +
  theme(axis.text.x  = element_text(angle=90)) +
  geom_hline(yintercept = 0, linetype = "dashed") + 
  labs(x="1 = Very liberal to 5 = Very Conservative", y = "Average Effect", caption = "", fill = "Ideology:") + 
  theme(legend.position = "none",axis.text.x = element_text(angle=37,size=12))

# Get rid of Other 
plot1$data %>%
  filter(value == "Boost2") %>%
  mutate(value = recode(value, "Boost2" = "Boost")) -> plot1$data


plot2 <- interplot(m = mod7, var1 = "treatment", var2 = "Ideology_num",facet_labs = c("Encourage"),
                   point = T,ci = .90) +
  steph_theme +
  theme(axis.text.x  = element_text(angle=90)) +
  geom_hline(yintercept = 0, linetype = "dashed") + 
  labs(x="1 = Very liberal to 5 = Very Conservative", y = "Average Effect", caption = "", fill = "Ideology:") + 
  theme(legend.position = "none",axis.text.x = element_text(angle=37,size=12))

# Get rid of Other 
plot2$data %>%
  filter(value == "Encourage2") %>%
  mutate(value = recode(value, "Encourage2" = "Encourage")) -> plot2$data


plot3 <- interplot(m = mod8, var1 = "treatment", var2 = "Ideology_num",facet_labs = c("Mandate Healthcare"),
                   point = T,ci = .90) +
  steph_theme +
  theme(axis.text.x  = element_text(angle=90)) +
  geom_hline(yintercept = 0, linetype = "dashed") + 
  labs(x="1 = Very liberal to 5 = Very Conservative", y = "Average Effect", caption = "TESS 2022. 90% Confidence Intervals", fill = "Ideology:") + 
  theme(legend.position = "none",axis.text.x = element_text(angle=37,size=12))

# Get rid of Other 
plot3$data %>%
  filter(value == "Mandate Healthcare2") %>%
  mutate(value = recode(value, "Mandate Healthcare2" = "Mandate Healthcare")) -> plot3$data


plot4 <- interplot(m = mod8, var1 = "treatment", var2 = "Ideology_num",facet_labs = c("Mandate Universities"),
                   point = T,ci = .90) +
  steph_theme +
  theme(axis.text.x  = element_text(angle=90)) +
  geom_hline(yintercept = 0, linetype = "dashed") + 
  labs(x="1 = Very liberal to 5 = Very Conservative", y = "Average Effect", caption = "", fill = "Ideology:") + 
  theme(legend.position = "none",axis.text.x = element_text(angle=37,size=12))

# Get rid of Other 
plot4$data %>%
  filter(value == "Mandate Universities2") %>%
  mutate(value = recode(value, "Mandate Universities2" = "Mandate Universities")) -> plot4$data

ggpubr::ggarrange(plot1,plot2,plot3,plot4)

# ---------------------------------------------------------------------------------------------------------------
# SI Fig. 2. Change in Vaccine Intentions and Attitudes Given Exposure to the Republican Voter Condition Relative to the Control Condition by Ideology (90% Confidence Intervals)
# ---------------------------------------------------------------------------------------------------------------

mod2 <- lm(Q4B_Likely_to_SecondDose100 ~ treatment * Ideology_num, data = df, weights = WEIGHT)
mod3 <- lm(Q4C_Likely_to_Boost100 ~ treatment * Ideology_num,  data = df, weights = WEIGHT)
mod8 <- lm(Q7_Support_Mandate_Caretakers ~ treatment * Ideology_num, data = df, weights = WEIGHT)
mod9 <- lm(Q8_Support_Mandate_Universities ~ treatment * Ideology_num, data = df, weights = WEIGHT)

plot1 <- interplot(m = mod2, var1 = "treatment", var2 = "Ideology_num",facet_labs = c("Second Dose"),
                   point = T,ci = .90) +
  steph_theme +
  theme(axis.text.x  = element_text(angle=90)) +
  geom_hline(yintercept = 0, linetype = "dashed") + 
  labs(x="1 = Very liberal to 5 = Very Conservative", y = "Average Effect", caption = "", fill = "Ideology:") + 
  theme(legend.position = "none",axis.text.x = element_text(angle=37,size=12))

# Get rid of Other 
plot1$data %>%
  filter(value == "Second Dose3") %>%
  mutate(value = recode(value, "Second Dose3" = "Second Dose")) -> plot1$data

plot2 <- interplot(m = mod3, var1 = "treatment", var2 = "Ideology_num",facet_labs = c("Boost"),
                   point = T,ci = .90) +
  steph_theme +
  theme(axis.text.x  = element_text(angle=90)) +
  geom_hline(yintercept = 0, linetype = "dashed") + 
  labs(x="1 = Very liberal to 5 = Very Conservative", y = "Average Effect", caption = "", fill = "Ideology:") + 
  theme(legend.position = "none",axis.text.x = element_text(angle=37,size=12))

# Get rid of Other 
plot2$data %>%
  filter(value == "Boost3") %>%
  mutate(value = recode(value, "Boost3" = "Boost")) -> plot2$data

plot3 <- interplot(m = mod8, var1 = "treatment", var2 = "Ideology_num",facet_labs = c("Mandate Healthcare"),
                   point = T,ci = .90) +
  steph_theme +
  theme(axis.text.x  = element_text(angle=90)) +
  geom_hline(yintercept = 0, linetype = "dashed") + 
  labs(x="1 = Very liberal to 5 = Very Conservative", y = "Average Effect", caption = "TESS 2022. 90% Confidence Intervals", fill = "Ideology:") + 
  theme(legend.position = "none",axis.text.x = element_text(angle=37,size=12))

# Get rid of Other 
plot3$data %>%
  filter(value == "Mandate Healthcare3") %>%
  mutate(value = recode(value, "Mandate Healthcare3" = "Mandate Healthcare")) -> plot3$data

plot4 <- interplot(m = mod8, var1 = "treatment", var2 = "Ideology_num",facet_labs = c("Mandate Universities"),
                   point = T,ci = .90) +
  steph_theme +
  theme(axis.text.x  = element_text(angle=90)) +
  geom_hline(yintercept = 0, linetype = "dashed") + 
  labs(x="1 = Very liberal to 5 = Very Conservative", y = "Average Effect", caption = "", fill = "Ideology:") + 
  theme(legend.position = "none",axis.text.x = element_text(angle=37,size=12))

# Get rid of Other 
plot4$data %>%
  filter(value == "Mandate Universities3") %>%
  mutate(value = recode(value, "Mandate Universities3" = "Mandate Universities")) -> plot4$data

ggpubr::ggarrange(plot1,plot2,plot3,plot4)

